clear

use "1850_1860_merged.dta"

set seed 349132

preserve
keep if age1850>=16 & age1850<=18 & male1850==1
keep if (white1850==1)

*Drop if dad's location not found
keep if occ1950_pop1850~=.
drop if realprop_parents1850==.

*Drop if child's occupation is missing, and drop if father's occupation is missing
drop if occ19501860==999
drop if occ19501860==997
drop if occ19501860==979

drop if occ1950_pop1850==999
drop if occ1950_pop1850==997
drop if occ1950_pop1850==979

tab age1850
tab age1860

egen taghistid = tag(histid1850)

sort age1850 taghistid
by age1850 taghistid: gen scale = _N if taghistid==1
by age1850 taghistid: egen sonrankx = rank(occscore1860) if taghistid==1
by age1850 taghistid: egen dadrankx = rank(occscore_pop1850) if taghistid==1
replace sonrankx = sonrankx/scale
replace dadrankx = dadrankx/scale

bysort histid1850: egen sonrank = max(sonrankx)
bysort histid1850: egen dadrank = max(dadrankx)

egen index = group(cty_fips1850)
su index
local indexmax = r(max)
gen constantreg = . 
gen coefficient = . 
forvalues i = 1(1)`indexmax'{
cap noisily reghdfe sonrank dadrank if index==`i' [pw=weight]

mat C = e(b)
replace constantreg = C[1,2] if index==`i' & _rc==0
su constantreg if index==`i'
replace coefficient = C[1,1] if index==`i' & _rc==0
su coefficient if index==`i'
}

gen p25 = constantreg + coefficient*.25
gen p50 = constantreg + coefficient*.5
gen p75 = constantreg + coefficient*.75

#delimit ;
local vars constant;
foreach x of local vars{;
gen `x'wt = `x'*weight;
drop `x';
rename `x'wt `x';
};

#delimit ;
collapse (mean) statefip1850 hasnormalschool p25 p50 p75 (sum) constant, by(cty_fips1850);


#delimit cr
gen constantrd = round(constant,1)
gen below10 = constantrd<10
count if below10==1
bysort statefip hasnormalschool: egen totctytype = count(hasnormalschool)
bysort statefip hasnormalschool: egen totbelow10 = total(below10)
gen fracbelow10 = totbelow10/totctytype
bysort statefip: egen maxfracbelow10 = max(fracbelow10)

tab statefip if maxfracbelow10>=.5
drop if maxfracbelow10>=.5

*Drop the states that had less than 10 in the specifications separating by parental assets
drop if statefip==12|statefip==19|statefip==48|statefip==55

tab statefip1850

#delimit cr
sort cty_fips

local vars p25 p50 p75
foreach x of local vars{

reghdfe `x' hasnormalschool, absorb(statefip1850) vce(robust)
su `x' if e(sample)==1 & hasnormalschool==0
local mean = r(mean)
tab statefip1850 if e(sample)==1
outreg2 using censustree1850mearank.xls, append excel dec(3) addstat(dvmean, `mean') adec(3)

#delimit ;
permute hasnormalschool _b[hasnormalschool], strata(statefip1850) reps(1000): reg `x' hasnormalschool i.statefip1850;
#delimit cr
}
